1 Téléchargement et préparation des données

1.1 Téléchargement des données

library(osmdata)
library(sf)
library(units)

# define a bounding box
q0 <- opq(bbox = c( 2.2247, 48.8188, 2.4611, 48.9019)) 

# extract Paris boundaries
q1 <- add_osm_feature(opq = q0, key = 'name', value = "Paris", value_exact = TRUE)
res1 <- osmdata_sf(q1)
paris <- st_geometry(res1$osm_multipolygons[1,])

# extract La Seine 
q2 <- add_osm_feature(opq = q0, key = 'name', value = "La Seine", value_exact = TRUE)
res2 <- osmdata_sf(q2)
seine1 <- st_geometry(res2$osm_multilines)
q2b <- add_osm_feature(opq = q0, key = 'name', value = "La Seine - Bras de la Monnaie", value_exact = FALSE)
res2b <- osmdata_sf(q2b)
seine2 <- st_geometry(res2b$osm_lines)

# extract Parcs and Cimetierres
q3 <- add_osm_feature(opq = q0, key = 'leisure', value = "park", value_exact = TRUE)
res3 <- osmdata_sf(q3)
parc1 <- st_geometry(res3$osm_polygons)
parc2 <- st_geometry(res3$osm_multipolygons)
q4 <- add_osm_feature(opq = q0, key = 'landuse', value = "cemetery", value_exact = TRUE)
res4 <- osmdata_sf(q4)
parc3 <- st_geometry(res4$osm_polygons)

# extract Quartiers
q5 <- add_osm_feature(opq = q0, key = 'admin_level', value = "10", value_exact = TRUE)
res5 <- osmdata_sf(q5)
quartier <- res5$osm_multipolygons

# extract Bars, Pubs & Cafe
q6 <- add_osm_feature(opq = q0, key = 'amenity', value = "bar", value_exact = TRUE)
res6 <- osmdata_sf(q6)
bar <- res6$osm_points
q7 <- add_osm_feature(opq = q0, key = 'amenity', value = "cafe", value_exact = TRUE)
res7 <- osmdata_sf(q7)
cafe <- res7$osm_points
q8 <- add_osm_feature(opq = q0, key = 'amenity', value = "pub", value_exact = TRUE)
res8 <- osmdata_sf(q8)
pub <- res8$osm_points

# extract Metro Stations
q9 <- add_osm_feature(opq = q0, key = 'station', value = "subway", value_exact = TRUE)
res9 <- osmdata_sf(q9)
metro <- res9$osm_points

# extract Restaurant
q10 <- add_osm_feature(opq = q0, key = 'amenity', value = "restaurant", value_exact = TRUE)
res10 <- osmdata_sf(q10)
resto <- res10$osm_points

1.2 Préparation des données

# use Lambert 93 projection
parc1 <- st_transform(parc1, 2154)
parc2 <- st_transform(parc2, 2154)
parc3 <- st_transform(parc3, 2154)
paris <- st_transform(paris, 2154)
seine1 <- st_transform(seine1, 2154)
seine2 <- st_transform(seine2, 2154)
quartier <- st_transform(quartier, 2154)
bar <- st_transform(bar, 2154)
pub <- st_transform(pub, 2154)
cafe <- st_transform(cafe, 2154)
metro <- st_transform(metro, 2154)
resto <- st_transform(resto, 2154)



# make layers pretty
parc <- do.call(c, list(parc1, parc2, parc3))
parc <- st_union(x = st_buffer(parc,0), by_feature = F)
parc <- st_cast(parc, "POLYGON")
parc <- parc[st_area(parc)>=set_units(10000, "m^2")]
parc <- st_intersection(x = parc, y = paris)

seine <- st_intersection(x = seine1, y = paris)
seine <- c(st_cast(seine[1])[2:5], seine[2])
seine <-c(seine, seine2)

bar <- bar[!is.na(bar$name),]
pub <- pub[!is.na(pub$name),]
cafe <- cafe[!is.na(cafe$name),]
bars <- c(st_geometry(bar), st_geometry(pub))
# bars <- c(st_geometry(bar), st_geometry(cafe), st_geometry(pub))
bars <- st_intersection(x = bars, y = paris)




quartier <- quartier[substr(quartier$ref.INSEE,1,2)==75,]

metro <- st_intersection(x = metro, y = paris)
metro <- metro[metro$station%in%"subway", ]
metro <- metro[metro$railway%in%"station",]
metro <- metro[metro$STIF.zone%in%1,]

resto <- resto[resto$amenity%in%"restaurant",]
resto <- resto[!is.na(resto$name),]
resto <- resto[,"cuisine"]
resto <- st_intersection(x = resto, y = paris)
resto$cuisine <- as.character(resto$cuisine)


save(list= c("paris", "quartier", "seine", "parc", "bars", "metro", "resto"), 
     file = "../data/paname.RData", compress = "xz")

2 Analyse et Cartographie des bars à Paris

2.1 Affichage des données

library(spatstat)
library(sf)
library(cartography)
library(maptools)
library(raster)
library(units)

load("../data/paname.RData")


# 1st map
par(mar = c(0,0,1.2,0))
plot(paris, col = "#D9D0C9", border = NA, bg = "#FBEDDA")
plot(parc, col = "#CDEBB2", border = NA, add=T)
plot(seine, col = "#AAD3DF", add=T, lwd = 4)
plot(st_geometry(quartier), col = NA,lty = 2, lwd = 0.2, add=T)
plot(bars, add=T, col = "#0000ff", pch = 20, cex = 0.2)
layoutLayer(title = "Cafés, bars et pubs à Paris", scale = 1,
            north = T,
            tabtitle = TRUE, 
            sources = "Map data © OpenStreetMap contributors, under CC BY SA.", 
            author = "cartography 2.0.2")

2.2 Nombre de bars par quartiers

# compute some data
quartier$nbar <- lengths(st_covers(quartier, bars))
quartier$dbar <- quartier$nbar / set_units(st_area(quartier), "km^2")


# 2nd map
par(mar = c(0,0,1.2,0))
plot(paris, col = "#D9D0C9", border = NA, bg = "#FBEDDA")
plot(parc, col = "#CDEBB2", border = NA, add=T)
plot(seine, col = "#AAD3DF", add=T, lwd = 4)
plot(st_geometry(quartier), col = NA,lty = 2, lwd = 0.2, add=T)
propSymbolsLayer(x = quartier, var = "nbar", inches = 0.2, col = "red", 
                 legend.pos = "topright", 
                 legend.title.txt = "Nombre de bars\npar quartier")
layoutLayer(title = "Cafés, bars et pubs à Paris", scale = 1,
            tabtitle = TRUE, 
            sources = "Map data © OpenStreetMap contributors, under CC BY SA.", 
            author = "cartography 2.0.2") 

# 3rd map
par(mar = c(0,0,1.2,0))
plot(paris, col = "#D9D0C9", border = NA, bg = "#FBEDDA")
choroLayer(quartier, var = "dbar", border = NA, 
           method = "quantile", nclass = 6,
           col = carto.pal("wine.pal", 6),
           legend.pos = "topright",
           legend.title.txt = "Densité de bars\npar km2", add = TRUE)
plot(parc, col = "#E2CCB5", border = NA, add=T)
plot(seine, col = "#AAD3DF", add=T, lwd = 4)
plot(st_geometry(quartier), col = NA,lty = 2, lwd = 0.2, add=T)
layoutLayer(title = "Cafés, bars et pubs à Paris", scale = 1,
            tabtitle = TRUE, 
            sources = "Map data © OpenStreetMap contributors, under CC BY SA.", 
            author = "cartography 2.0.2") 

# 2nd map
par(mar = c(0,0,1.2,0))
plot(paris, col = "#D9D0C9", border = NA, bg = "#FBEDDA")
plot(parc, col = "#CDEBB2", border = NA, add=T)
plot(seine, col = "#AAD3DF", add=T, lwd = 4)
plot(st_geometry(quartier), col = NA,lty = 2, lwd = 0.2, add=T)
propSymbolsChoroLayer(x = quartier, var = "nbar", var2 = "dbar",
                      inches = 0.2,
                      method = "quantile", nclass = 6,
                      col = carto.pal("wine.pal", 6),
                 legend.var2.pos  = "topright", 
                 legend.var2.title.txt = "Densité de bars\npar km2",
                 legend.var.title.txt = "Nombre de bars\npar quartier")
layoutLayer(title = "Cafés, bars et pubs à Paris", scale = 1,
            tabtitle = TRUE, 
            sources = "Map data © OpenStreetMap contributors, under CC BY SA.", 
            author = "cartography 2.0.2")

2.3 Analyse lissée de la répartition des bars

bb <- as(bars, "Spatial")
bbowin <- as.owin(as(paris, "Spatial"))
pts <- coordinates(bb)
p <- ppp(pts[,1], pts[,2], window=bbowin)
ds <- density.ppp(p, sigma = 200, eps = c(20,20))
rasdens <- raster(ds) * 1000 * 1000
rasdens <- rasdens+1
par(mar = c(0,0,1.2,0))
bks <- getBreaks(values(rasdens), nclass = 12, method = "equal")
cols <- colorRampPalette(c("black","#940000", "white"))(length(bks)-1)
plot(paris, col = NA, border = NA, main="", bg = "#FBEDDA")
plot(rasdens, breaks= bks, col=cols, add = T,legend=F)
legendChoro(pos = "topright",cex = 0.7,
            title.txt = "Densité de Bars\ndans un voisinage gaussien\n(bars par km2)",
            breaks = bks-1, nodata = FALSE,
            col = cols)
plot(seine, col = "#AAD3DF", add=T, lwd = 4)
plot(parc, col = "#CDEBB235", border = NA, add=T)
plot(quartier$geometry, add=T, border = "white",lty = 2, lwd = 0.1)
plot(bars, add=T, col = "#0000ff90", pch = 20, cex = 0.1)
plot(paris, col = NA, add=T)
layoutLayer(title = "Cafés, bars et pubs à Paris", scale = 1,
            tabtitle = TRUE, 
            sources = "Map data © OpenStreetMap contributors, under CC BY SA.", 
            author = "cartography 2.0.2")

3 Analyse et cartographie des restaurants à Paris

3.1 Les données

# 1st map
par(mar = c(0,0,1.2,0))
plot(paris, col = "#D9D0C9", border = NA, bg = "#FBEDDA")
plot(parc, col = "#CDEBB2", border = NA, add=T)
plot(seine, col = "#AAD3DF", add=T, lwd = 4)
plot(st_geometry(quartier), col = NA,lty = 2, lwd = 0.2, add=T)
plot(resto, add=T, col = "blue", pch = 20, cex = 0.1)
layoutLayer(title = "Tous les restaurants", scale = 1,
            north = T,
            tabtitle = TRUE, 
            sources = "Map data © OpenStreetMap contributors, under CC BY SA.", 
            author = "cartography 2.0.2") 

# nf <- layout(mat = matrix(data = c(1,1,2,1,1,3,4,5,6), nrow = 3, byrow = T), 
#              widths = c(1,1,1), respect = FALSE)

3.2 Cartes de densités

densityMap <- function(x = resto, cuisine = NULL, title, sigma ){
  opar <- par(mar = c(0,0,0,0))
  if(!is.null(cuisine)){
    x <- x[x$cuisine %in% cuisine, ]
  }
  bb <- as(x, "Spatial")
  bbowin <- as.owin(as(paris, "Spatial"))
  pts <- coordinates(bb)
  p <- ppp(pts[,1], pts[,2], window=bbowin)
  ds <- density.ppp(p, sigma = sigma, eps = c(20,20))
  rasdens <- raster(ds) * 1000 * 1000
  rasdens <- rasdens+1
  xx <- getBreaks(values(rasdens), nclass = 12, method = "equal")
  cols <- colorRampPalette(c("black", "#940000", "white"))(length(xx)-1)
  plot(paris, col = NA, border = NA, main="")
  image(rasdens, breaks= xx, col=cols, add = T,legend=F)
  legendChoro(pos = "topright",cex = 0.7,title.cex = 1,
              title.txt = paste0(title, "\nn=",nrow(pts)),
              breaks = xx-1, nodata = FALSE,values.rnd = 0,
              col = cols)
  plot(seine, col = "#AAD3DF", add=T, lwd = 4)
  plot(parc, col = "#CDEBB235", border = NA, add=T)
  plot(quartier$geometry, add=T, border = "white",lty = 2, lwd = 0.1)
  plot(st_geometry(x), add=T, col = "blue", pch = 20, cex = 0.1)
  barscale(size = 1)
  mtext(text = "cartography 2.0.2\nMap data © OpenStreetMap contributors, under CC BY SA.", 
        side = 1, line = -1, adj = c(0.01,0.01), cex = 0.6, font = 3)
  north(pos = c(661171.8,6858051))
  par(opar)
}

densityMap(x = resto, title = "Tous les restaurants", sigma = 500)

densityMap(x = resto, cuisine = c('japanese','sushi'), title = "Retaurant japonais", sigma = 500)

densityMap(x = resto, cuisine = c('pizza','italian'), title = "Restaurant italiens", sigma = 500)

densityMap(x = resto, cuisine = c('chinese'), title = "Retaurant chinois", sigma = 500)

densityMap(x = resto, cuisine = c('korean'), title = "Restaurant coréens", sigma = 500)

densityMap(x = resto, cuisine = c('indian'), title = "Restaurant indiens", sigma = 500)

densityMap(x = resto, cuisine = c('asian', "thai", "vietnamese"), title = "Restaurant asiatiques", sigma = 500)

# 06
# detach("package:spatstat",character.only=TRUE,unload=TRUE)
library(sf)
x <- st_distance(metro, metro)
diag(x) <- 10000
dmin <- apply(x, 2, min)
metro$nbar <- lengths(st_covers(st_buffer(x = metro, dist = 375), bars))
png("Paris06.png", width = 1200, height = 663, res = 100)
par(mar = c(0,0,1.2,0))
plot(paris, col = "#D9D0C9", border = NA, bg = "#FBEDDA")
plot(parc, col = "#CDEBB2", border = NA, add=T)
plot(seine, col = "#AAD3DF", add=T, lwd = 4)
plot(st_geometry(quartier), col = NA,lty = 2, lwd = 0.2, add=T)
x <- st_buffer(x = metro, dist = 375)
x$r <- dmin/2
choroLayer(x = x, var = "nbar", col = carto.pal("wine.pal",5),
                       method = "fisher-jenks", nclass = 5, add=T,
                 legend.pos = "topright", 
                 legend.title.txt = "Nombre de bars à proximité\nde la station de métro")
labelLayer(x = metro[x$nbar>50,], txt = "name", overlap = F, 
           show.lines = F, halo = T, r = 0.09)
layoutLayer(title = "Cafés, bars et pubs à Paris", scale = 1,
            tabtitle = TRUE, 
            sources = "Map data © OpenStreetMap contributors, under CC BY SA.", 
            author = "cartography 2.0.2") 
dev.off()

# 7
library(osrm)
bars <- st_sf(bars)
options(osrm.server = "http://0.0.0.0:5001/", osrm.profile = "walk")

distA2 <- osrmTable(src = as(metro, "Spatial"),
                    dst = as(bars, "Spatial"))
mat <- distA2$durations

metro$nb10 <- apply(mat,1, FUN = function(x){sum(x<=5)})


par(mar = c(0,0,1.2,0))
plot(paris, col = "#D9D0C9", border = NA, bg = "#FBEDDA")
plot(parc, col = "#CDEBB2", border = NA, add=T)
plot(seine, col = "#AAD3DF", add=T, lwd = 4)
plot(st_geometry(quartier), col = NA,lty = 2, lwd = 0.2, add=T)
metrob <- st_buffer(metro, 150)
choroLayer(metrob, var = "nb10",  add=T, legend.pos = "topright",
           method = "equal", nclass = 7, col = carto.pal("red.pal", 7),
           legend.title.txt = "nb bars à moins\nde 5 min à pied",
           legend.values.rnd = -1,
           border = "grey20", lwd = .5)
layoutLayer(title = "Cafés, bars et pubs à Paris", scale = 1,theme = 'wine.pal',
            tabtitle = TRUE, 
            sources = "Map data © OpenStreetMap contributors, under CC BY SA.", 
            author = "cartography 2.0.2")